Heterogeneous shear in hard sphere glasses 
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There is growing evidence that the flow of driven amorphous solids is not homogeneous, even if the 
macroscopic stress is constant across the system. Via event driven molecular dynamics simulations 
of a hard sphere glass, we provide the first direct evidence for a correlation between the fluctuations 
of the local volume-fraction and the fluctuations of the local shear rate. Higher shear rates do 
preferentially occur at regions of lower density and vice versa. The temporal behavior of fluctuations 
is governed by a characteristic time scale, which, when measured in units of strain, is independent 
of shear rate in the investigated range. Interestingly, the correlation volume is also roughly constant 
for the same range of shear rates. A possible connection between these two observations is discussed. 



PACS numbers: 

Introduction. — Heterogeneous flow and shear band- 
ing are central to the rheology of complex fluids and are 
widely observed in many industrial and natural materi- 
als, such as foams, emulsions, pastes, or even rocks (HQ- 
Despite its ubiquitous appearance, many aspects of this 
phenomenon are still not well-understood. In the sim- 
plest case, shear banding can be captured by a nonmono- 
tonic dependence of the shear stress, a, on the shear rate, 
7 For certain complex fluids, like colloidal gels, 

shear banding can be associated to the competition be- 
tween a structural phase transition and a shear [B[ . How- 
ever, in other systems, such as dense hard sphere (HS) 
colloidal suspensions and granular materials 0, flow 
heterogeneity is often observed without such accompany- 
ing structural changes. The rhcological response of these 
systems is essentially determined by the competition be- 
tween an inherent slow dynamics and the acceleration 
caused by the external drive || Q . This may lead to a 
spatially and temporally heterogeneous flow if the system 
is close to the yielding threshold 

Recently, it has been proposed that flow localization in 
dense hard sphere suspensions can be alternatively ratio- 
nalized in terms of shear-concentration coupling (SCC) 
Q, a well-known feedback mechanism for flow instabil- 
ity in complex fluids fl3| . In this model, regions of high 
(low) shear rate are associated with high (low) diffusiv- 
ity, giving rise to a shear-induced flux from high to low 
shear-rate regions. This leads to an increase (decrease) 
of density and thus viscosity in low (high) shear-rate re- 
gions, thereby enhancing shear-rate fluctuations further. 
This, in turn, enhances the migration of particles from 
high to low shear-rate regions, unless balanced by diffu- 
sive counter flux. As pointed out recently enhance- 
ment of density fluctuations and the associated formation 
of voids can also be essential for the occurrence of frac- 
ture. 

However, while the experimental data in Q are inter- 
preted consistently within the proposed macroscopic pic- 



ture, no test of the basic underlying assumptions, such 
as the presence of a correlation between shear-rate and 
concentration fluctuations or the growth of emergent ve- 
locity fluctuations, has been provided so far for colloidal 
hard sphere glasses. This is not surprising, as the rele- 
vant density fluctuations that trigger the initial instabil- 
ity are quite small and hardly accessible to experiments. 
Furthermore, the available experimental time window for 
the observation of velocity fluctuations is limited to a few 
hundred percent strain, thus making a temporal analysis 
rather difficult. 

Here, we study these and related issues via event- 
driven molecular dynamics (MD) simulations of a poly- 
disperse hard sphere system. It is explicitly shown that 
fluctuations of the local volume fraction, 8(f), are corre- 
lated to the fluctuations of the local shear rate 5j. More 
precisely, a decrease of local density is accompanied by 
an increase of the local shear rate and vice versa. How- 
ever, while the relative amplitude of shear-rate fluctua- 
tions increases upon decreasing shear rate, the magni- 
tude of the fluctuations of volume fraction remain con- 
stant. The temporal behavior of the fluctuations is char- 
acterized by a dominant time scale which is identical for 
both observables 6<p and S'j. When measured in units 
of strain, this time scale is of the order of a few hun- 
dred percent strain and fairly independent of the im- 
posed shear rate in the studied range. Interestingly, the 
correlation volume, as determined from the maximum 
of the four-point susceptibility of density fluctuations, 
Xa = ^\{fq{t)) ~ {fq(t)) 2 ]i is a lso roughly constant in 
the investigated range (N is the particle number and 
f q (t) = N^ 1 Yli=i exp[iq • (r.;(t) - n(0))] is the incoher- 
ent scattering function at wave vector q). This suggests a 
possible link between an inhomogeneous shear and a dy- 
namically heterogeneous environment with a given spa- 
tial correlation. 

Simulation setup. — We perform event driven MD 
simulations of a polydisperse (11%) hard sphere system 
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in 3D [15J. The volume fraction is varied from below to 
above the glass transition point, which, for the present 
polydisperse svjstem, is located at a volume fraction of 
<j) g ps 58.5% The quiescent properties of this sys- 

tem have been studied extensively in (T3). The temper- 
ature is fixed at T = 1 via velocity rescaling. In all 
the simulations reported here, N = 7000 particles [28| 
are placed in a random configuration between two walls. 
The walls are made of the same kind of particles as the 
bulk, but have infinite mass. A constant shear is applied 
by moving the walls with velocities ±J7 wa ii in the ±x di- 
rection (planar Couette flow). The shear gradient points 
in the +z direction. The (overall) shear rate is defined 
as 7 a v = 2U wa \i/Lz, where Lz is the distance between 
the walls. Local quantities, such as velocity and density 
profiles, are computed as an average over particles within 
distinct layers of finite width parallel to the walls. Note 
that, for simplicity of notation, we will drop the sub- 
script "av." Local quantities will be identified by their 
arguments, e.g., j(z) and (ft(z). 

Theoretical description. — Central to the theory of 
shear-concentration coupling are a non-Newtonian con- 
stitutive relation for the shear stress and the notion of 
a noncquilibrium particle pressure [l3j |. Both quanti- 
ties are coupled through the Navier-Stokes equations. In 
Figs. Ufa) andQJb), we show the dependence of the shear 
stress on shear rate and volume fraction in the yield stress 
regime, as obtained from our simulations. In the glassy 
phase ((ft > 0.585 <-> $ > 0.873), our data can be well 
described by a Hcrschcl-Bulkley expression 
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where p, n, o~q and A are constants to be determined by a 
fit. $ = 4>/4> m is the reduced volume fraction ((ft m = 0.67 
corresponds to random close packing). The first term in 
Eq. (JTJ represents the dynamic yield stress while 
the second term accounts for the effect of shear. The dy- 
namic yield stress can be naturally associated with shear 
heterogeneity, if the rigid regions are understood to be 
locally below the yieldin g th reshold, while the liquidlikc 
regions are above 0, Isl. TloL Il9j |. However, it does not 
explain the mechanism of the initial flow instability. 

The particle pressure obtained from our simulations for 
different densities and shear rates is shown in Figs. QJc) 
and [ljd). Noting that the pressure data show similar 
behavior to the shear stress, it appears natural to write 
(m, r, IIo and B are fit parameters) 
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The first term constitutes the pressure contribution re- 
sponsible for ordinary concentration diffusion (20| . The 
shear-rate dependence of the particle pressure is a mani- 
festation of shear-induced particle migration, also known 
as "dilatancy" (HHH. 



Within SCC, the feedback responsible for heteroge- 
neous flow occurs if an initial density excess leads, via 
d^a, to a locally lower shear rate, which in turn drives, 
via c^n, further particle migration towards the low shear 
region. The flow instability develops as soon as these 
fluxes cannot be overcome anymore by the concentration 
and viscous momentum diffusion, described by c^II and 



djcr, respectively. A detailed calculation 
that this happens if 



f = (^n)(a*(r)/(5*n)(^<7) > 1. 



shows 
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In the interesting lower shear-rate regime, F reduces to 
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As seen from Figs. Ufa) andQJc), f° r the glassy phase 
in the limit of low shear rates, shear stress as well as 
pressure show quite a similar dependence on shear rate. 
Indeed, trying various fit procedures revealed that m ss n 
leads to consistent fit results for all the simulated data. 
This suggests that Fq is practically independent on 7. 
Using this information in F(j c , <J>) = 1, the critical shear 
rate for the onset of instability is obtained: 
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Heterogeneous shear. — Figure O^a) shows the veloc- 
ity profile at successive strain intervals, illustrating the 
occurrence of flow heterogeneity in the studied polydis- 
perse HS model in the glassy phase (0 = 0.6075). Note 
the qualitative similarity to experiments of colloidal hard 
spheres 6] and to simulations of a binary Lenard-Jones 
glass [10(. Here, we go a step further and perform a de- 
tailed survey of temporal fluctuations of both local shear 
rate, Sj — 7(2) — 7, and volume fraction, 5(j) = <fi(z) — (ft 
[Fig. Hlb)]. The data suggest that a positive 8 (ft is often 
accompanied by a negative ^7 and vice versa. Moreover, 
even though the temporal analysis is performed at con- 
siderably higher 7 (due to limited computer time), the 
time scale of fluctuations in Figs.^a) and[^b) seems to 
be roughly the same in units of strain. 

We compute the statistical average of the instanta- 
neous correlations between fluctuations of the local shear 
rate and volume fraction, C < j 1 ^{z) = (8 , y(z)8(ft(z) / ((fry)) . 
A typical result on Cj,p(z) is shown in Fig. [3Ja) for an 
average shear rate of 7 = 10~ 5 . We first note that most 
data points are distributed around a constant value in- 
dicated by a dashed line in the plot. This value is small 
but definitely nonzero (see the error bars, representing a 
standard deviation determined from independent runs). 
It is worth mentioning that the computational effort to 
obtain the present statistical accuracy is quite significant 
(40 independent runs, each with a duration of 60 days on 
a 3GHz CPU). The large values of C^{z) at z ps ±9 are 
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FIG. 1: Shear stress a and particle pressure II versus (a),(c) 
imposed shear rate 7 and (b),(d) reduced volume fraction 
$. Solid lines are fits to (a),(b) Eq. (JTJ) with parameters 
p ~ 2.355, n = 0.4, a = 0.01 19, and A = 18 - 35 and 
to (c),(d) Eq. Q with parameters r = 3.8 — 4.1, m = 0.4, 
n ~ 8.4 - 10.05, and B = 0.0015 - 0.007. From bottom to 
top, different curves correspond to (a),(c) $ = 0.8060(<^> = 
0.5400), 0.8657(0.5800), 0.8970(0.6014), 0.9067(0.6075), 
0.9189(0.6157), 0.9251(0.6198) and 0.9299(0.6230) and (b),(d) 
7 = 10" 5 ,4 x 10" 5 , 10" 4 ,4 x 10" 4 , 10" 3 ,4 x 10" 3 and 10" 2 . 

a consequence of wall effects. Indeed, the volume fraction 
close to the walls is slightly (by about 1%) but systemat- 
ically larger than the average volume fraction, leading to 
a corresponding decrease of wall shear rate. A detailed 
study of this issue will be presented elsewhere. For the 
reminder of this Letter, we will be concerned with bulk 
properties only. 

To test the stability phase diagram, Eq. Eq. ([5]), via 
our simulations, we identify a heterogeneous flow by re- 
quiring max|(i>(z) — jz)/U wa a\ > 0.18. As can be seen 
from Fig. |3jb), simulation and theoretical predictions are 
in reasonable agreement. Thus, SCC seems to describe 
at least the onset of instability quite well for the present 




FIG. 3: (a) Correlation between fluctuations of local shear 
rate and volume fraction, determined within parallel layers 
of two-particle diameter thickness, (b) Comparing simulation 
results to Eq. ([5]) with p = 2.355, r = 4.1, m, n = 0.4, A = 
34.5 and B = 0.007. The dashed line marks the parameter 
range for which a detailed temporal analysis of fluctuations is 
performed (see Fig. [4}. 

system. 

Time scale of fluctuations. — To elucidate the nature 
of the instability, we study layer-resolved fluctuations, 
8 A f(z,t) and 5<f>(z,t), and determine, via a Fourier anal- 
ysis, the peak frequency and the associated character- 
istic time scale of fluctuations at <f> = 0.6075 for shear 
rates ranging from homogeneous to inhomogeneous flow 
regimes [dashed line in Fig. Efb)] . We first note that this 
time scale is independent of z if the distance to the walls 
is above a few particles in diameter (data not shown). 
More importantly, it is practically identical for both the 
observables £7 and 5(f> [Fig. Ufa)]. Moreover, when mea- 
sured in the units of strain, it is independent of applied 
shear in the investigated range. Using the same set of 
data, we also determine the relative amplitude of fluctu- 
ations as a function of shear rate [Fig. Eta)]. In marked 
contrast to the behavior of the characteristic time scale, 
this quantity depends on the specific observable. In par- 
ticular, Sj/j increases significantly with decreasing im- 
posed shear rate, while 6<fi/(j> is perfectly constant in the 
interesting limit of low shear rates, i.e., when entering 
the inhomogeneous flow regime. 

The above observations have a number of consequences 
for a possible interpretation of flow heterogeneity. First, 
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FIG. 2: (a) Fluctuation of the velocity profile in a sheared 
HS glass. Each curve is computed as an average over the 
indicated strain interval (100% strain corresponds to a time 
of I/7). (b) <57/7 and 54>/4> versus time at the center of the 
simulation cell for a 300 times higher shear rate. 
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FIG. 4: (a) Dependence on shear rate of the relative ampli- 
tude and the characteristic time scale of fluctuations, (b) The 
function X4 whose maximum is a measure of the correlation 
volume. 
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the presence of a correlation between fluctuations of local 
shear rate and volume fraction [FigJHJa)] and the analy- 
sis of the stability diagram [Fig. EJb)] are, at least in the 
first sight, in favor of SCC theory. However, the basic 
assumption of SCC theory is that fluctuations of shear 
rate give rise to fluctuations of volume fraction and vice 
versa. This implies that, as the instability is approached, 
both Sj and Scj) should grow to some extent. This is 
at odds with the data shown in Fig. 0Ja). To see this, 
we first note that, in contrast to the supercooled state, 
where a slight change of volume fraction can lead to con- 
siderable variations of viscosity as the glass transition is 
approached, the viscosity of a glass changes only little 
with volume fraction. This can easily be inferred from 
a comparison of shear stress for volume fractions below 
and above 4> g = 0.585 ($ g = 0.873) in Fig. (Ha). It is thus 
not plausible that tiny (in the sense of hardly detectable) 
changes of volume fraction can account for the observed 
increase of shear-rate fluctuation amplitude. 

In this light, it is tempting to search for a link be- 
tween the observed temporal behavior of Sj and S(j) and 
dynamic heterogeneity. Without any claim for rigor, we 
start from the intuitive idea that some regions in the ma- 
terial are more mobile (i.e., appear more fluidlike) than 
others and therefore can support larger shear. These re- 
gions are not static. Rather, they continuously form and 
dissociate. The time scale of the fluctuations in the flow 
velocity is thus expected to reflect a time scale inherent 
to this microscopic dynamics. This inherent time scale 
should in turn be related to the characteristic size of these 
dynamically correlated regions. If this idea is consistent, 
the constant time scale [29!] observed in our simulations 
would then imply a constant correlation volume. This is 
indeed borne out in Fig. |U(b), where we plot X4> whose 
maximum is a measure of the correlation volume [24| , as 
a function of shear rate. We remark that dynamic corre- 
lations are not visible in the fluctuations of one-particle 
density so that, based on this idea, no conclusion can be 
made as to how local density fluctuations should behave. 
Therefore, the different behaviors of the fluctuation am- 
plitudes of shear rate and density shown in Fig. 0k are 
not a priori in conflict with this picture. 

In conclusion, we find that, although the theory of 
shear-concentration coupling describes some aspects of 
flow heterogeneity in hard sphere glasses, serious incon- 
sistencies remain, such as a significant increase of the 
relative amplitude of shear-rate fluctuations upon ap- 
proaching the instability, as opposed to constant fluc- 
tuations of volume fraction. A more stringent test of the 
SCC theory would be to numerically solve the underlying 
Navier-Stokes equations using the ex per imental (or sim- 
ulated) cr(7,(/)) and 11(7,0) as input [25] . An important 
question to be answered would be whether the original 
SCC theory predicts stable shear bands or a fluctuat- 
ing behavior as observed in our simulations. Further- 
more, keeping in mind that solvent effects on the system 



response are negligible in the limit of high <f> and low 
7, it would also be interesting to work out the conse- 
quences of shear-concentration coupling by solving the 
compressible Navier-Stokes equation (complemented by 
a non-Newtonian constitutive law) only for the colloids. 
When seeking for alternative routes, on the other hand, 
it must be kept in mind that, in its present form, the 
proposed SCC theory for a HS glass is a hydrodynamic 
approach with a fully local constitutive law. Dynamic 
correlations, however, give rise to nonlocal effects. These 
could possibly be taken into account (e.g., by means of a 
correlation length of some measure of "fluidity" 26[) by 
introducing generalized hydrodynamic equations. 
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